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TRANSFORMATION OF TWO AND THREE-DIMENSIONAL 
REGIONS BY ELLIPTIC SYSTEMS 

A major portion of our effort during this period has been in transferring 
our computational work from the LRC computer to the IRIS Graphics Workstation 
at MSU. Our first successful application was in the computation of a 
conservative solution of a simple hyperbolic equation on an overlapping grid. 
That example is included in the attached report. The report is a revision of 
work included in our last status report, and the results will be presented at 
the First International Conference on Numerical Grid Generation. As our 
experience on the IRIS increases, more complicated geometric configurations 
will be considered. Several conclusions concerning computations on 
overlapping grids are apparent. Problems only occur when there is a major 
difference in grid spacing on the individual component grids. In the case of 
hyperbolic equations, it is necessary that both interpolation and 
extrapolation be applied at the grid boundaries. When interpolated values are 
used at outflow boundary points, excessive oscillations in the numerical 
solution may be the result. The same conclusions would be valid for more 
complicated systems of hyperbolic equations such as the Euler equations for 
inviscid flow. Some of the solution values would be extrapolated at the 
overlap boundary, the exact number depending on the number of characteristics 
pointing out of the overlap region. It is also possible that similar boundary 
conditions may be needed for some parabolic equations such as high Reynolds 
number viscous flow equations. 

Considerable effort has been expended on the development of 
three-dimensional conservative interpolation procedures. While the 
theoretical development is a straight forward extension of the two-dimensional 
concepts, the technical difficulties in implementing a feasible algorithm 



appear to be overwhelming. Consequently, our work on three-dimensional 
problems will be limited to the case where the overlap boundaries coincide 
with grid surfaces. Further progress in three-dimensions will await the final 
algorithm development and verification of the general two-dimensional method. 

We have begun our investigation of grid smoothing procedures during this 
reporting period. Although no significant new results can be reported, the 
direction of future research has been established. It has been decided that 
the first grid smoothing algorithms will be based on the concepts of 
variational grid generation. In recent years there have appeared several 
modifications of the basic variational method, but in all cases the 
fundamental idea is to control geometric grid quantities such as distance 
between points, angle of intersection of grid lines, and cell volumes. Thus 
consideration is given to all grid properties effecting error in finite 
difference or finite volume computations. All of the variational methods give 
rise to systems of nonlinear equations which are lengthy and often difficult 
to solve. It is this property that we wish to avoid by applying the 
variational principle on a local basis. The algorithms will be simpler 
because the objectives are limited. The objective is not to generate a 
completely new grid but to improve an existing algebraic grid. The grid 
smoothing methods will be applied in the same spirit as the smoothing filters 
commonly used in data analysis. The smoothing algorithm will only be applied 
a few times with little interest in the eventual convergence of the grid or 
properties of a converged grid. 



INTERFACE PROCEDURES FOR OVERLAPPING GRIDS* 


C. Wayne Mastin 
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Abstract 


Interpolation at grid boundaries is studied for the 
purpose of solving partial differential equations using 
either implicit or conservative explicit finite-difference 
methods on multi -component overlapping grid systems. 


1 . INTRODUCTION A multi-component grid system, in which 
several computation grids are used, is required in the nu- 
merical solution of many fluid dynamics problems involving 
flow within or about a complicated geometric configuration. 
From a grid construction point of view, the simplest proce- 
dure is to generate each component grid independently with a 
sufficient overlap so that information can be transmitted 
from one grid to the other. The development and analysis of 
solution procedures on this type of grid system was studied 
by Starius [8,9]. The practical application of the method to 
the solution of problems in computational fluid dynamics was 
demonstrated in the papers by Atta [1], Atta and Vadyak [2], 
and Thompson [11]. This was followed by further studies on 
interpolation techniques by Kreiss [5] and Mastin and 
McConnaughey [6]. Each successful application, such as the 
recent results of Steger and Buning [10] and Benek, et al 
[3]. serves to reinforce the need for additional work on the 
implementation of numerical methods on overlapping grids. 


*Research supported by NASA Langley Research Center 
under Grant No. NSG 1577. 



ORIGINAL PAGE IS 
OF POOR QUALITY 


Two popular numerical methods In computational fluid 
dynamics are the alternating direction implicit algorithms 
and the explicit algorithms derived for the solution of 
conservation laws. Both methods have been used in the solu- 
tion of problems on composite grid systems. However, in each 
case there are properties of the numerical solution which are 
lost when information is transmitted between the individual 
grids. With the implicit algorithms, there is no technique 
for generating advanced solution values at all boundary 
points of each component grid. Lagging some of the boundary 
point values can lead to a loss of accuracy in the solution 
of transient problems. Of course, the temporal step length 
could be reduced, but that would defeat the purpose of choos- 
ing an Implicit method. It is also possible that lagging may 
effect the stability of the method, although no problems of 
that kind have been reported. 

A conservative finite-difference scheme is often se- 
lected when solving partial differential equations in conser- 
vation form. Since the classical interpolation formulas were 
not derived with conservation properties in mind, their use 
in finite-difference approximations on composite grids would 
result in the I033 of an exact conservation property. An 
interpolation scheme for conservative finite-difference 
methods was first proposed by Rai [ 7 ]. He considered compos- 
ite grid systems which did not overlap but Joined along 
common grid lines. Berger [M] indicated how the method of 
Rai could be generalized and extended to overlapping grids. 

This report will describe ways of eliminating the time 
lag in implicit solutions and will present a general algo- 
rithm for constructing conservative interface conditions. 

The variables x and y are used as the spatial variables in 
the partial differential equations. Since, in the general 
case, one would need to solve the equations on a curvilinear 
grid, all equations would be transformed to curvilinear 
coordinates before applying the finite-difference algorithm. 
In each development the partial differential equation is 
sufficiently generally so that the method can be applied to 
the original or the transformed equation without change. For 
simplicity, mixed derivative terms which are normally lagged 
in the single grid case and source terms are not included in 
the partial differential equations. 


2 . IMPLICIT METHODS The fundamental concepts are quite 
simple and can be demonstrated by considering the one-dimen- 
sional equation 


Uj. - Lu 


( 1 ) 


where L is a differential operator of the form 



Lu - Au x + Bu xx 


Since implicit methods generally require linearization of the 
difference equations, it may as well be assumed that L is 
linear. Suppose that two grids G., and G 2 are given on the 
intervals [a,d] and [c,b], respectively, where 


a < c < d < b. 


If M denotes the usual second-order difference approximation 
of L, then the Crank-Nicolson equation can be written as 


u 


n+1 

i 


u 


n 



<MuJ +1 


♦ MUj*) . 


( 2 ) 


Here i is the spatial index, n is the temporal index, and At 
is the step length. Now suppose values on G^ and G 2 are 
known at level n and values at level n+1 are to be computed 
on G 1 . While solution values needed in (2) at x-d can be 
interpolated from G 2 for level n, the corresponding values at 
level n+1 are unavailable. If these unknown values are 
replaced by the values at level n, then the local truncation 
error at the neighboring interior point is Increased by a 
term on the order of 0(At 1 2 3 /Ax 2 ). The value Ax represents the 
spatial grid spacing on Gj, or the spacing at x-d In the case 
of a nonuniform grid. In any event, when Ax is small, this 
lagging of solution values will seriously degrade the tempo- 
ral accuracy of the approximation. The error can be reduced 
by following a particular order in updating the solution 
values at the interior grid boundary points. The correct 
sequence of computations is indicated in the following steps. 


1 . Calculate u n+ ] 

2. Calculate u n+ ; 

3. Calculate u n+2 
M. Calculate u n+2 


on G 1 with level 
on G 2 with level 
on G 2 with level 
on Gj with level 


n values at x-d. 
n+1 values at x-c. 
n+1 values at x-c. 
n+2 values at x-d. 


Now the error induced by using the previous value at x-d in 
step 1 is offset by the use of the advanced value in step 4. 
In fact, the local truncation error at the neighboring inte- 
rior point is increased by a term of order 0(At’/Ax 2 ) when 
the solution is advanced from level n to n+2. The same error 
reduction would also occur at x-c. 



Clearly, this four-step alternating grid scheme is only 
a partial solution. Unless the solution exhibited a linear 
growth or decay, there would still be points with a local 
truncation error of order one whenever At-Ax. However, this 
does not necessarily mean that the global error in the nu- 
merical solution would be increased to that order. The 
actual error in the solution would also depend on other 
factors such as the extent of the overlap. Note that the 
same updating procedure could be applied to implicit methods 
other than the Crank-Nicholson method, but the reduction in 
local truncation error would not be the same. 

The alternating grid concept has also been used in the 
development of another method for implementing implicit 
algorithms on composite grid systems. This method also alter- 
nately employs the forward difference explicit equation 


uj +1 - u£ ♦ At MuJ (3) 

and the backward difference implicit equation 

uj +1 - u} ♦ At Mu^ +1 • (4) 


The computational sequence is illustrated in the following 
four-step procedure which advances the solution from level n 
to level n+2. 

1. Calculate u n+ | on Gj using (3). 

2. Calculate u n+ ^ on G 2 using (4). 

3. Calculate u n+2 on G 2 using (3). 

4. Calculate u n+2 on G^ using (4). 


The method alternates the explicit and implicit calculations 
in the same manner as the well-known hopscotch algorithm. 
Thus, the name hopscotch will be associated with this method. 
The method has several desireable properties. All values 
needed at the grid boundaries c and d can be computed by 
interpolation from solution values at the correct time level. 
The overall method is second-order accurate in time and 
unconditionally stable. This fact follows by noting that the 
combined sequence of (3) followed by (4) is equivalent to a 
Crank-Nicolson step with step length of 2At. 

The consequences of lagging solution values at grid 
boundaries can be even more serious for multi-dimensional 
problems. Suppose for example that the operator L in (1) is 



defined as 


Lu - Au x + BUy ♦ Cu xx 4 Du yy . 

Let M x and M y denote the difference approximations of the x 
and y derivative parts of L. If the parabolic equation (1) 
is solved by an ADI method, such as Peaceman-Rachford, the 
algorithm becomes 


u n*1 /2 . u n . |£. (H]i u n *' /2 . M y u”) . (5a) 

u"*' . u n+1/2 . |i (M]tU n*1/2 „ HyU n *1). (5b) 


Now the error that occurs in the first step of the algorithm 
is further magnified in the second step. This argument can 
be made more precise by noting that the Peaceman-Rachford ADI 
method is a perturbation of the two-dimensional Crank- 
Nicolson method with a perturbation term 


M^(u n - u n+1 ) . 


A lagged value in this term produces a truncation error term 
on the order of 0(At^/Ax 2 Ay 2 ) . The alternating grid procedure 
would reduce this to 0(At l */Ax 2 Ay 2 ). The one-dimensional 
hopscotch algorithm would not be a computationally efficient 
method for two-dimensional problems. However, the same 
effect can. be realized by inserting additional steps in the 
ADI algorithm. The procedure is again demonstrated using two 
grids and G 2 . Note that equation (5a) can be written as 


y n+1/2 _ u n + + M y u n ) (6a) 

u n+1 /2 . yh+1/2 4 ( „ xU n4l/2 _ ^n^ (6b) 


while (5b) can be replaced by 


V n+ 1 - u n+1/2 ♦ (V"* 172 ♦ M y u n+1/2 ) 

u n+ 1 - V n+1 ♦ (MyU n+1 - M y u n+1/2 ) . 


(6c) 

(6d) 


These split forms would require additional computations, and 
should only be used to generate interpolated values at inte- 
rior grid boundaries. The following steps illustrate one 
possible method of computation. 


1. Calculate v n+ | ^ on Gj using (6a) 

2. Calculate u n+ ^ 2 on G 2 using (5a) 

3. Calculate u n+ |' 2 on Gj using (6b) 
M. Calculate v I1+ ' on G 2 using (6c) 

5. Calculate u n+1 on G^ using (5b) 

6. Calculate u n+1 on G 2 using (6b) 


Note that at each step the necessary boundary values for one 
grid can be Interpolated from values at the correct level on 
the other grid. 

The efficacy of the alternating grid and hopscotch 
methods is exhibited in the solution of a one-dimensional 
model problem. The parabolic equation 

u t ♦ (u-c)u x - p u xx (7) 

has an exact solution 


u(x,t) 


1 

2 


(1 - tanh 


2(x+1 )+(2c-1 ) t 
‘•U 


). 


This equation is solved on the interval [-2,2] with the exact 
initial value at t-0 and boundary values at x--2 and x-2. A 
second-order linearization and the usual central difference 
approximations are used. An overlapping set of two grids on 
the Intervals [-2,. 125] and [-.125,2] is constructed. The 
solution for values of c-0.4 and y-0.05 is computed using 
three different methods. The form of the actual solution 
indicates that an increase in t would result in a translation 
of the graph in the positive x direction. When a numerical 
solution is computed with the Crank-Nicolson equation (2) and 
the values at x-±.125 are Jagged, there is a marked deviation 
between the numerical and analytic solutions as they pass 
through the overlap interval. Although the numerical solu- 
tion lags behind the actual solution, they are qualitatively 



similar with no indication of Instability in the numerical 
solution. A comparison of the solutions at various times is 
plotted in Figure 1 . The lag in the numerical solution is 
eliminated when the alternating grid method is used. A 
careful examination of Figure 2 reveals an anomaly in the 
graph at the grid points adjacent to the interior boundary 
points x-±.125. This is more evident on the enlargement in 
Figure 3. Note that the problem occurs only at the points 
where the exceptional difference approximation is employed. 
The most accurate numerical solution for this example is 
calculated using the hopscotch algorithm. That solution 
appears in Figure H. In all of these figures, linear inter- 
polation was used to determine solution values at grid bound- 
aries. 


3. METHODS FOR CONSERVATION LAWS If the conservation 
equation 


u t + [f(u)] x - 0 


( 8 ) 


is solved on a composite grid system, then there must be some 
means of transferring the flux f(u) from one grid to the 
other. There are two feasible alternatives. Either the 
solution u can be calculated by interpolation and then f(u) 
evaluated, or f(u) can be interpolated directly from the flux 
values on the other grid. The conservative difference sch- 
emes which will be discussed require interpolation of fluxes. 
However, before proceeding in that direction, a comparison of 
the two interpolation techniques will be included. 

Suppose a solution value u* at a boundary point of grid 
G 1 is computed by linear interpolation from the solution 
values u 1 _ 1 and Uj defined on grid G 2 . Then an interpolation 
formula of the form 


u* - au 1 _ 1 + BUj, o+B-1, 

holds, and the flux can be evaluated as f(u*). Now if Uq is 
the actual value of the solution at the boundary point of G. , 
and hence the true flux value is f(Urt), then the interpolati- 
on procedure introduces an error as is seen in the following 
expansion. 


f(ou 1 _ 1 +6u i ) - f(u 0 )+f u (u 0 )(au 1 _ 1 +3u i -u 0 ) 


+ . 









Since an expansion at the boundary point Xq yields 


auj-i * euj - u 0 “ 7 u xx 0 (a( Xi-rx 0 ) + 6(x r x 0 ) 2 ), 


+ ... 


the leading term in the local truncation error is 


7 f u (u 0 )u xx0 (o(x i-r x 0 )2 + BCxj-Xq) 2 ). 


Whenever the option of interpolating the fluxes is selected, 
then the boundary flux f* is calculated directly as 

f* - a ftu^) + 3f(u i ). 

Expanding about the solution Uq, and noting the additional 
second order term, 

o f(u i _ 1 ) ♦ B f(uj) - f ( U q ) + f u (u 0 )(ctu 1 _ 1 +Bu 1 -u 0 ) 

+ J f uu (u 0 )(a(u i-1~ u 0 )2 + 6(ui-u 0 ) 2 ) 

The leading term in the local truncation error now has the 
form 

2 ^ ^"u ^ **0 ^ ^xxO + ^uu^**0^ u xo^ a ^*i— 1 "**0^ + ®^i“Xq) )). 

It is clear that both procedures give an Interpolation error 
which is 0 ( Ax 2 ) . 


There are many conservative finite-difference algorithms 
for solving conservation laws of the type (8). Most of the 
basic algorithms of practical Interest can be written as 



g (u; +1 ,u}) - g (u{,uj. 1 ) 


(9) 


uj +1 “ U 8 - 


For notational convenience, let 

g i+1/2 “ g ^ u i+1 * u i^ * 

with the implication here being that 8^/2 is an approxima- 
tion at x 1+1 / 2 * Thia notation is appropriate for the central 
difference approximations such as the Lax or Lax-Wendroff 
schemes. When one-sided or upwind differencing is used, the 
fractional index 1+1/2 would be replaced by i or 1+1. 

Given a grid with grldpolnts x^, 1-0,1 1, the discrete 

conservation property states that 

1-1 1-1 

Y H n+1 . Y Ii n + !T n - <r n 
i u i i u i + g I-1/2 g 1/2 * 

i-1 i-1 

The same result can be obtained from (8) by using numerical 
integration from x 1 / 2 to Xt-i/ 2 871(1 the flux approximation 
determined by g. It is this derivation that will be used in 
the composite grid approach. 

Let and Go be grids defined on the Intervals [a,d] 
and [c,b], a<c<d<b. The grid G 1 has points x lt 1-0,1,..., I 
and grid spacing Ax, and G 2 has points y., J-0,1,...,J and 
spacing Ay. The difference equations will be written in 
terms of scaled solution values v and w defined by 


v - uAx and w - uAy 


On G^ , the difference equation has the form 

v n+1 _ v n _ h n _ h n 
V 1 v i h i+1/2 h i-1/2 


and on G 2 , 

w^ +1 " wj - k j+i/2 " k J-1/2 • 

There is good reason for writing the equations in this form. 
First of all, the grid spacing need not be Included in the 
interpolation formulas, but more importantly, this is the 
required form of the difference equations when computing on 
moving grids. 



The correct interface conditions can now be derived by 
extending the grid functions to piecewise linear functions 
and integrating. Suppose a value k 1 / 2 is needed. Then the 
interval [a,b] is partitioned into two subintervals Ca.yi/p] 
and Cy^/ 2 »b]. If y -\/2 lies in the interval C x i-i /2* x l+l /2^’ 
and hi_i /2 and h i+i/2 are known » then the value for k 1 /2 can 
be calculated from the integral property 


* 1/2 

h x 

x l/2 


y J-l/2 

k x “ k J-1/2 ” h 1/2 * 

y l/2 


Assuming that h and k are piecewise linear, it is easily seen 
that the needed value is the linear lnterpolant defined as 


k 1/2 " ° h l-1/2 + e h i+1/2 » 

where 

a « x iti i2iiM2 f g _ y.i 12. - # 

x i+1/2“ x i-1/2 x i+1/2 _x i-1/2 


By the same argument, the interval [a,b] can be partitioned 
into [a,Xj_.j/ 2 ] and C x j-i/ 2 » b ^ and the interpolation formula 
for the value hj_j/ 2 is 

h I-1/2 " a k J-1/2 + B k J+1/2 

where 

a _ y Jt17?~ x T-,1/2 ? g . x T-17?~yj-1/? . 

y J+1/2~yj-1/2 ’ yj+1/2 _y j-1/2 


The same linear Interpolation would be used on a non- 
uniform grid. Of course, the scaling factor would vary from 
point to point. An interpolation formula could also have 
been derived using the original equation (9), however the 
difference in grid spacing on and G 2 would have resulted 
in the appearance of a scaling factor In the interpolation 
formula. 

A modification of this approach can be used to develop 
conservative interface conditions for two and three-dimen- 
sional problems. The general two-dimensional conservation 
law is 



I 


u t + f x + 8y - °* 

where x and y are now the spatial variables. A difference 
approximation has the form 


,n+1 _ n 
r i.J v i.J 


♦ h 1 


n 


1+1/2.J 


h n * t,n _ i,n 

n i-1/2,j + k i,j+1/2 k i,j-1/2* 


( 10 ) 


The grid function v is the product of the solution u and the 
Jacobian (or cell area), and the values of h and k are, up to 
a scalar factor, flux values in the direction of the curvi- 
linear coordinate lines. Let and G 2 be overlapping grids 
and suppose values of h are required along the i - 1/2 grid 
line of . For now, it is assumed that the endpoints of the 
grid line are on the boundary of the physical region. The 
points of G 1 along the grid line are labeled using parameter 
values from any convenient parameterization. Thus, let 


P J " (x 1/2,J» *1/2, J>» J-0,1....,J, 

while points of intersection of the 1-1/2 grid line of G 1 
with all grid lines of G 2 are ordered and labeled 

*1 m (x i,J+fi’ or < x i*«,J* *i + 6.J ) ' 

1-0 ,1 , . . . ,L 


where 6 denotes a fractional index between 0 and 1, and i,J 
are the indices of some point in G 2 . The first step in the 
transfer of flux values from G 2 to G 1 is to define flux 
values at the points q^. If q^ lies on an 1-constant grid 
line, as in the first case above, then a value hi is computed 
by Interpolating the grid function k. On the other hand, if 
Qfc lies on a J -constant grid line, then h£ is computed by 
Interpolating h. Now the i - 1/2 grid line divides the 
physical region into two parts, one covered by G 1 and the 
other covered by a subset of G 2 . If piecewise linear flux 
functions are constructed along the grid lines in each sub- 
region and an Integration of the flux derivatives over the 
complete region is performed, then the conservation property 
requires that 




_ l . — 



2 


J-l 



+ 


M -* 


L-l 


i 

l»l 



( 11 ) 




Here the 1 and n indices in (10) have been suppressed. 

The Interpolation formula will be defined using a set of 
basis functions. Two additional parameter values are intro- 
duced by extrapolating from the parametric interval. Let q_ 1 
- 2q Q - and q L+ < - 2q L - q^-i • Let ^ be the piecewise 
linear function, with knots q^, l--1,0 # ...,L+1 f defined as 




! 1, m-i, 

0 , mtl 


where t-0,1,...,L, and m— 1 ,0,. . . t L+1 . The following inte- 
grals can be easily computed from the parametric values of 
the points along the grid line. 




q L+1 

+1 

q -1 


A t*0 " 


p 1/2 

*1 

JPo 


A 


l.J “ 


Pj+1 /2 
Pj-1 /2 


J-1,2 J-1 


A 




Pj 

h 

Pj-1/2 


These Integrals are used in calculating the coefficients of 
the interpolation formulas. The formulas can now be written 
as 


h|* for J-0 and J-J, 


y 


L A 


h 4 - 2 




l.J 


and 


J 0 V ht ' rorj ■’• 2 J '’ 


( 12 ) 


The fact that property (11) holds is readily verified. 




* $ n J ■ Jo Jo ‘V h * 

■ Jo h t7. J. 


i j-0 
h* L-l 

• * ih h i 


+ — 


A few remarks are sufficient to indicate how the same 
interpolation method can be employed in more general compos- 
ite grid configurations. If interpolation is required on a 
boundary component consisting of several i -constant and 
J -constant segments, then each boundary segment could be 
^ treated separately with either an h value or a k value calcu- 

lated from equations (12). However, the extrapolated parame- 
ter values would not be used in computing the coefficients, 
but Instead a single parameterization would be defined for 
the entire boundary component. If the boundary component 
were a clofeed contour in the interior of the physical region, 
then the special boundary Interpolation formulas for j-0 and 
J-J in (12) would be unnecessary. 

The selection of a set of piecewise linear basis func- 
tions to define the interpolation coefficients may be changed 
with only slight modification. One could just as easily use 
piecewise constant functions or use higher degree polynomials 
such as quadratics or cubics. The degree of interpolation 
may have differing effects on the numerical solution. The 
use of a piecewise constant basis may produce shock-like 
discontinuities, whereas a linear basis would tend to smear 
out any discontinuities in the solution. 

The conservative finite-difference scheme of MacCormack 
^ is used to solve the parabolic equation (7) which can be 

written in conservation form as 



This equation is solved on two overlapping grids on the 
intervals [-2,. 25] and [-.25,2] with the same values of c-O.M 
and y-0.05 as in the previous section. The numerical dissi- 
pation in the MacCormack scheme permits the use of a coarser 
grid than was used for the implicit methods. Figure 5 con- 
tains the solution plotted for various values of t. For this 
example, there was no noticeable difference between solutions 
computed with flux values interpolated and those computed 
with interpolation of solution values. 

While the choice of interpolation methods will have some 
effect on the numerical solution, a more fundamental question 
for hyperbolic equations is when interpolation should be used 
to generate boundary values for a grid and when the boundary 
values should be determined by interior solution values using 
some numerical boundary condition. If the characteristic 
direction at a boundary point of an overlap region is exterior 
to the region, then any attempt to impose boundary values by 
interpolation would result in an overdetermined problem. 

Since the same problem is being solved on both grids in the 
overlap region, the difference between the interpolated value 
and the value of the solution determined by the 
characteristics would be significant only for problems with 
shocks or other high gradient regions. When the difference 
is significant, as in the following example, the numerical 
solution exhibits the familiar oscillatory form. 



Figure 5. Conservative explicit solution 




A simple two-dimensional conservation law is given by 

p t + (up) x ♦ ( vp)y - 0. 

Two grids are used to construct a solution on a rectangular 
region. Let G 1 denote a grid covering the points with 
coordinates satisfying 

0 S x S 1, 0 S y S 1 , 


and let G 2 cover the points with 


0.9 S x S 1.9, OS y S 1. 


Now choose constant values u-1 and v«0. The necessary 
boundary values are given as 

! 1 .0 if yS0.5 , 

0.5 if y>0.5 . 

The steady-state solution to this essentially one-dimensional 
problem is p-1.0 for yS0.5 and p«0.5 for y>0.5. A numerical 
solution is computed using MacCormack's method with 
conservative interpolation formulas at the boundary x-1 of 
grid G^ and the boundary x-0.9 of a coarser grid G 2 . A plot 
of the solution at the interior grid points appears in Figure 
6. The error caused by the Inconsistency in the Interpolated 
value near the discontinuity is clearly evident. This 
particular solution is computed with an interpolation formula 
derived from piecewise linear basis functions. The same 
behavior is observed when piecewise constant basis functions 
are selected. In this example the error is practically 
eliminated by extrapolating from G^ to obtain boundary 
values of p along x-1 . The successful solution of the 
problem is illustrated in Figure 7. 

It is noted that even when numerical boundary conditions 
are used, the numerical solution still satisfies a 
conservation property provided the region can be partitioned 
along contours where the conservative Interpolation scheme is 
applied. In the above example the divergence Integral would 
be approximated using values of p on G^ for 0SxS0.9 and 
values on G 2 for 0.9SxS1. 9. 


A. CONCLUSIONS The accuracy of the transient solution of 
a hyperbolic or parabolic partial differential equation is 
dependent upon the procedures used to transfer information 


Figure 6. z=p(x,y) surface with interpolation on 
overlap boundary 



Figure 7. z*p(x,y) surface with Interpolation and 
extrapolation on overlap boundary 


between grids in a composite grid system. The error in the 
numerical solution can be reduced by using the techniques 
developed here. While the attempt has been to construct 
algorithms that are easy to implement, the degree of diffi- 
culty would ultimately be linked to the complexity of the 
grid structure. 
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